function [F] = CDKM(K1, F,K2,sita)
    fK1f = sum(F .* (K1 * F))';
    fK2f = sum(F .* (K2 * F))';
    ff = sum(F)';
    m_all = vec2ind(F')';
    obj(1) = sum(diag(sita{1}).*(fK1f ./ ff))+sum(diag(sita{2}).*(fK2f ./ ff));

    for iter = 2:50
          for i = 1:size(F, 1)
            m = m_all(i);
            if ff(m) == 1
                % avoid generating empty cluster
                continue;
            end

            Y_A = F' * K1(:, i);

            fK1f_s = fK1f + 2 * Y_A + K1(i, i);
            fK1f_s(m) = fK1f(m);
            ff_k1 = ff + 1;
            ff_k1(m) = ff(m);

            fK1f_0 = fK1f;
            fK1f_0(m) = fK1f(m) - 2 * Y_A(m) + K1(i, i);
            ffK1_0 = ff;
            ffK1_0(m) = ff(m) - 1;


            Y_B = F' * K2(:, i);

            fK2f_s = fK2f + 2 * Y_B + K2(i, i);
            fK2f_s(m) = fK2f(m);
            ff_k2 = ff_k1;
            ff_k2(m) = ff_k1(m);

            fK2f_0 = fK2f;
            fK2f_0(m) = fK2f(m) - 2 * Y_B(m) + K2(i, i);
            ffK2_0 = ff;
            ffK2_0(m) = ffK1_0(m);


            delta = diag(sita{1}^(-1)).*(fK1f_s ./ ff_k1 - fK1f_0 ./ ffK1_0)  + diag(sita{2}^(-1)).*(fK2f_s ./ ff_k2 - fK2f_0 ./ ffK2_0);

            [~, p] = max(delta);
            if p ~= m
                fK1f([m, p]) = [fK1f_0(m), fK1f_s(p)];
                ff([m, p]) = [ffK1_0(m), ff_k1(p)];

                fK2f([m, p]) = [fK2f_0(m), fK2f_s(p)];
                %ff([m, p]) = [ffK2_0(m), ff_k2(p)];

                F(i, [p, m]) = [1, 0];
                m_all(i) = p;
            end
          end
        obj(iter) = sum(diag(sita{1}).*(fK1f ./ ff))+sum(diag(sita{2}).*(fK2f ./ ff));

        if iter > 2 && (obj(iter) - obj(iter - 1)) / obj(iter - 1) < 1e-9
            break;
        end
        
end
